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ABSTRACT 

Recently much work has been done with the goal of establishing ultra-low emission 
aircraft gas turbine combustors. A significant portion of this effort is the development of 
three-dimensional computational combustor models. The KIVA-II computer code which 
is based on the ICEd-ALE numerical scheme is one of the codes selected by NASA to 
achieve these goals. 

This report involves a simulation of jet injection through slanted slots within the 
Rich-burn/Quick-quench/Lean-burn (RQL) baseline experimental rig. The RQL combustor 
distinguishes three regions of combustion. This work specifically focuses on modeling the 
quick-quench mixer region in which secondary injection air is introduced radially through 
12 equally spaced slots around the mixer circumference. 

Steady state solutions are achieved with modifications to the KIVA-II program. 
Work currently underway will evaluate thermal mixing as a function of injection air velocity 
and angle of inclination of the slots. 
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NOMENCLATURE 


DR* = jet-to-mainstream density ratio. 

= Tm/Tj 

i,j,k = grid coordinates 

i - plane coordinates ranging between 
the centerline and outer radius, 
in the radial direction, 
j - plane coordinates ranging between 
the azimuthal boundaries of 
the sector mesh, in the azimuthal 
direction respectively, 
k - plane coordinates ranging from the 
base of the mesh to the outlet, in 
the axial direction. 

j* = jet-to-mainstream momentum flux ratio. 

= [ (Vj) (Vj) (RHOj) ] / [ (Vm) (Vm) (RHOm) ] 

NX = the number of cells in the radial 

dircetion, with NXP = NX + l. 

NY = the number of cells in the azimuthal 

direction, with NYP = NY + 1. 

r* = jet-to-mainstream velocity ratio. 

= (Vj/Vm) 

RHOj *» jet density (g/cubic cm) . 

RHOm = mainstream density (g/cubic cm) . 

Tj = jet temperature (°K). 

Tm = mainstream temperature (°K). 

U, V, W = cartesian velocity components (cm/sec) , 

W denoting the axial component. 

Vj = jet inflow velocity (cm/sec) . 

Vm = mainstream inflow velocity (cm/sec) . 

Win = axial velocity applied at k = 1 plane. 

X,Y,Z = cartesian coordinates, with Z 

indicating downstream distance from 
k = 1 plane. 
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CHAPTER I 


INTRODUCTION 

I. 1. Motivation 

Much of the research currently underway with gas 
turbine combustion is in the development of ultra-low 
emission aircraft gas turbine combustors, to ensure no 
significant ozone depletion due to emissions from future 
civil transport aircraft. Several designs are being 
investigated to meet this goal. One design is the Lean/ 
Premixed-Prevaporized (LPP) combustor. This combustor avoids 
droplet combustion and operates with a low fuel/air ratio in 
the primary zone to achieve low flame temperature with 
resulting low NOx emissions. Another gas turbine combustor 
design that shows promise is the Rich-burn/ Quick-quench/ 
Lean-burn (RQL) combustor. The RQL baseline combustor 
utilizes slanted slot injection orifices in the quench 
region. Research with an experimental prototype of the RQL 
combustor is currently underway at the NASA Lewis Research 
Center in Cleveland, Ohio. The baseline combustor attempts 
to separate combustion into three distinct burn zones (rich, 
intermediate, lean) . In the quench region, which is the 
intermediate combustion zone, air is injected into the 
combustion products as they leave the rich burn region. The 
injected air is at ambient temperature and is directed 
radially inward from slanted slots located along the 
cylindrical combustor wall. 
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With the development of civil transport aircraft 
traveling at Mach 3. 0-5.0 stringent emissions controls are 
becoming increasingly important. The new civil transport 
aircraft will cruise high in the stratosphere where 
emissions may easily interact with ozone. Ozone depletion 
occurs by the following reaction with nitrogen oxide [ 1 ]: 

(1) NO + 0 3 N0 3 + 0 3 

(2) NO 2 +0 •+■ NO + 0 2 

Equation (1) shows destruction of ozone ( O 3 ) and the 
second equation represents reconstitution of nitrogen oxide. 
Lefebvre [1] breaks down nitrogen oxide sources into three 
different categories: 

1. Thermal NO 

2 . Prompt NO 

3 . NO contained within the fuel 

Nitrogen oxide is generated in high temperature regions of 
the combustor and is most prominent above 1800° Kelvin. 
Prompt NO occurs directly within the flame front and is a 
result of interaction among intermediate species during the 
combustion process. Fuel NO is generated from organic 
nitrogen contained within its molecular structure. Thermal 
NO is the source of NO that both the LPP Combustor and the 
rqL Combustor are designed to minimize by maintaining well 
controlled low flame temperatures during the combustion 
process. Thermal NO is a direct result of the extended 
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Zeldovich chain mechanism expressed by the following 
equations : 


(1) 

0 

+ 

n 2 

N + NO 

(2) 

0 

+ 

NO 

«-*■ 0 2 + N 

(3) 

N 

+ 

OH 

«-* H + NO 


The RQL Combustor maintains low flame temperature with a 
high equivalence ratio (ER) in the rich zone. The combustor 
has a design ER of 1.6. The RQL combustor also reduces the 
production of NO from nitrogen in the fuel since higher 
equivalence ratios and lower flame temperatures both 
minimize this conversion rate. 

Combustion products leaving the rich zone of the RQL 
Combustor enter a converging region separating the rich and 
mixing regions. The mixing region is equipped with twelve 
injection orifices that inject the secondary combustion air 
to complete the combustion process. An experimentally 
verified, bulk swirl velocity is predicted to develop 
downstream of the slanted slots. The combustion process is 
enhanced in this region by both an additional supply of 
oxygen and highly turbulent mixing that exists due to the 
jet-mainstream interaction. The mixture travels through a 
radial expansion as it flows into the lean zone. Figure 
(1-1) shows a schematic of the RQL experimental facility. 
Many choices exist for the size and shape of the injection 
orifices in the mixing region, as well as the choice of a 
single or double row of orifices. 
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1.2. Literature Survey 

The mixing region of the combustor is in essence an 
example of turbulent jet mixing in a confined crossflow. 

Many papers are available that describe penetration of a 
circular jet into an open crossflow. However, only a few 
discuss jet penetration into a confined crossflow. Opposed 
jet injection into a confined crossflow has several specific 
applications in currently operating gas turbine combustors, 
such as air injected into the dilution zone to shape outlet 
temperature profiles from the combustor. Another example of 
this type of flow exists with some combustors in the primary 
zone immediately downstream of the injector. Air is injected 
normal to the axial flow inducing a recirculation zone. The 
recirculation zone entrains unburned combustion products 
aiding in the control of the combustion reaction rate and 
soot formation [1]. The RQL's mixing section introduces a 
new application for opposed jet injection within the gas 
turbine combustor. 

Previous work describing the velocity and temperature 
profiles of a circular jet into an open cross stream has 
been done by Keffer & Baines [2] and Ramsey and Goldstein 
[3]. More work describing near field measurements of the jet 
trajectory, temperature and velocity profiles was conducted 
by Moussa et al. [4]. This work was limited to a circular 
jet. Prediction of the distortion of a circular jet into a 
kidney shaped profile was made possible by Adler & Baron [5] 
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with an involved application of integral methods. A finite 
difference scheme utilizing a two equation turbulence model 
was used by Patankar et al. [6] to verify the experiments 
of Kefer & Baines and Ramsey & Goldstein. 

Experiments conducted by Kamotani & Greber [7] produced 
velocity and temperature measurements for a single jet 
injected into a confined cross stream. Jones & McGuirk used 
a finite difference scheme with a two equation turbulence 
model to reproduce the experimental observations of Kamotani 
& Greber. 

Holdeman has conducted considerable theoretical and 
experimental work to investigate the confined crossflow jet 
injection problem. Holdeman has made several important 
observations of jet mixing within a confined crossflow 
[8,9]. These are some of the important conclusions of his 
work that pertain directly to the RQL Combustor: 

1 . Mixing improves with increasing 
downstream distance from the jets. 

2. Momentum flux ratio (J*) is the most 
significant flow field variable. 

3. Variation of the density ratio at a 
specified momentum flux ratio yields 
minimal variation of the flow field. 

4. Smaller momentum flux ratios require 
a greater downstream distance for 
equivalent mixing. 

5 . Flow area convergence improves the 
quality of mixing. 
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Holdeman ' s numerical work was done with an interactive 
computer code written in Basic, which was developed to aid 
in modeling dilution zone air injection, which is used to 
shape the outlet temperature profile of gas turbine 
combustors [9]. His experimental work includes work with 
many orifice geometries, including 45 degree slots of the 
type found in the RQL baseline geometry. Holdeman found 45 
degree slots, in a rectangular duct, produce lower mixing 
quality and less penetration into the mainstream of a 
confined cross-flow compared to equivalent area streamlined 
slots or circular orifices [10]. 

1.3. Current Study 

The objective of the current research is to determine 
orifice geometries and operating conditions for optimum 
mixing in the quick-quench region. Experimental work done by 
the UTRC with the RQL prototype, shown in figure (1-1), 
evaluated three momentum flux ratios at a constant density 
ratio with three different 45 degree slot geometries having 
varying slot length and width. This work utilized planar 
digital imaging with an oil aerosol marking the jet. The 
study concluded optimum mixing occurred at the intermediate 
momentum flux ratio ( J* ) of 17.8. The flowfield in the 
vicinity of the mixing region of the RQL combustor is very 
complex. Vortices are found to exist between the jets near 
the slot inlet. For parallel slots, two vortices of equal 
size but counter rotating exist between the slots. However, 
for a 45 degree slanted slot one clockwise rotating vortex 
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is dominant in size, while the other is virtually 
nonexistent. As jet velocity increases, radial jet 
penetration and axial penetration of the swirling flow into 
the lean zone both increase. The swirling profile may be a 
significant contributor to the quality of mixing if optimum 
mixing in the experimental apparatus occurred at J* = 17.8 
and not the higher momentum flux ratio of 78.0. This 
conclusion is based upon the dispersion of oil droplets and 
not thermal mixing. 

An analysis of the jet-to-mainstream interaction for 
this mixer requires a three-dimensional numerical model. 
KIVA-II with its well established ICEd-ALE finite volume 
numerical scheme is employed to investigate the flow field. 
KIVA-II uses the standard k- e turbulence model. The thermal 
mixing performance will be determined from the predicted 
temperature and velocity fields for specified orifice 
geometry and operating conditions. 


7 




Figure 1-1 
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CHAPTER II 


PROGRAM DEVELOPEMENT 
II. 1. Program Description 

All rapid mixer geometry nonreactive flow mixing 
calculations were conducted with KIVA-II, a reactive flow 
code developed at Los Alamos National Laboratory. This code 
is capable of solving transient, two- and three-dimensional, 
chemically reactive fluid flows with sprays. KIVA-II solves 
the unsteady equations of motion for a turbulent, chemically 
reactive mixture of ideal gases, coupled with equations for 
a single component vaporizing fuel spray. The code resolves 
spatial differences by subdividing the computational region 
into small hexahedrons. Coordinates of the hexahedron cell 
corners (vertices) are calculated as functions of time, 
which allows for a Lagrangian, Eulerian, or combined 
description of the flow. A major advantage of the KIVA-II 
code is its ability to use nonorthogonal cells. 

Nonorthogonal cells are required in the converging and 
diverging regions of the baseline geometry. The ICE method, 
which is used by KIVA— II, utilizes an implicit solution for 
pressure. This allows larger time steps than purely explicit 
differencing of the pressure gradient terms. 

Each time step of the solution may involve both a 
Lagrangian and a Eulerian rezone phase. Vertices move with 
the fluid in the Lagrangian phase, and then are repositioned 
to user specified positions during the rezone phase. After 
repositioning the vertices, the flowfield is interpolated 
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onto the new computational mesh. Explicit methods are used 
to calculate convection in the rezone phase. 

KIVA-II uses the k- e turbulence model modified to 
include volumetric expansion effects and spray-turbulence 
interactions. Boundary layer drag and wall heat transfer are 
calculated with a modified turbulent law-of-the-wall model. 
II. 2. Governing Equations 

KIVA-II is capable of solving both laminar and 
turbulent flows. For turbulent flows transport coefficients 
are derived from turbulent diffusivity. The governing 
equations solved by the code are as follows: 

Mass Conservation Equation 


3P. 


•c 

m 


_E + $. (p S) = + p‘ 

Vh m m i 

for the mth species, where P m is the mass density of 

species m, u is the fluid velocity, is the diffusive 

mass flux of species m, and p^ is the rate of change of 

due to chemical reaction. The diffusive mass flux is given 

by Fick's Law: 

J = pD^Y 
m m 

where p is the total mass density and Y m is the species 
mass fraction. The mass diffusivity D , which is assumed to 
be the same for all species, is determined using the Schmidt 
number Sc, (assumed constant) : 

D = 


pSc 


where y is viscosity. 
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1^- + (pu) = 0 


Momentum Cohservation Equation 

- + $-(puu) = -^p - ^ ^ pk + V-o 

where p is the fluid pressure, and k is the specific 
turbulent kinetic energy. The viscous stress tensor is 
Newtonian and is expressed as: 

c = y [ Vu + (Vu) T ] + XV *uU 


where P and A are the shear and bulk viscosity 
coefficients, U is the unit dyadic tensor, and 
superscript T denotes the transpose. The viscosity 
coefficients are given by: 


p 


air 


+ 


c p 
P 



£ 


and 


A = 


2 

'J y 


where e is the specific turbulent kinetic energy 
dissipation rate and is an empirical turbulence 
constant equalling 0.09. A Sutherland formula is used 
for the molecular viscosity of air with and A 2 

constant: 


A 1 T 


3/2 


air 


T + A, 
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Energy Conservation Equation 


^§ 2 . + ?.<p 3 l) - -p$.u - V-q + pe + Q c 


where I is the specific internal energy and Q C is a 
source term due to chemical heat release. The total heat 
flux q , is gained by summing the effects of heat conduction 
and enthalpy transport by mass diffusion: 


- 5 C + 5 d = -K?T - |h m J : 


m 


where K is the thermal conductivity, T is the fluid 
temperature, and is the species specific enthalpy. The 

thermal conductivity is determined using the Prandtl number 
Pr (assumed constant) . 


K = 



where is the specific heat at constant pressure. 
Turbulence Equations 


^ • (puk) + j pk^*u = a :^u + vk] - pe 

r k 

(pue) + [i c + c ]pe^*u = 

0 L ^ e i C 3 

$e] + £ [c a :^u - c pe] 
e K e l ~ e 2 


where c £ ^ , c £ ^ , c £ ^ , Pr k , and Pr £ are empirical 

constants having the following values: 


c =1.44 

c = 1.92 

e l 

e 2 

Pr k - 1.0 

Pr = 1.3 
e 



- 1.0 
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Equations of state assume ideal gas relationships. The 
following are the equations of state used by KIVA-II. 


p = R TE( p /W ) 
r o m m 

I (T) =E( p /p )I (T) 
m m 

c (T) =Z( p/p )c (T) 
p m pm 

h (T) = I (T) + R T/W 
m m o m 


W 

m 

I (T) 
m 

c _ C T ) 

pm 

h (T) 
m 


m 


the universal gas constant 

molecular weight of species m 

specific internal energy of species m 

specific heat at constant pressure of 
species m 

enthalpy of species m 
density of species m 


II. 3. Mesh Generation 

KIVA-II was originally set up to allow for inflow along 
the bottom (k=l) boundary, and outflow along the right 
(i=nxp) and top (k=nzp) boundaries. To incorporate a slanted 
slot for secondary air injection into the mixing section the 
code was modified to allow inflow along a portion of the 
right boundary. To model the 45 degree slanted slot, k and j 
planes were spaced so looking at the right face of the 
cylindrical mesh along the i = nxp plane a line connecting 
the vertices could be drawn 45 degrees from the horizontal. 
Figure (2-la) shows a schematic description of an actual 45 
degree slot, its dimensions and important program variables 
used in generation of the mesh within a preprocessor. Figure 
(2-lb) shows dimensions of the baseline mixer from k=l 
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to 60 . 


For ns equally spaced slots, the flowfield will be 
periodic in 2( u ) /ns degree intervals. Therefore, an 
important capability of KIVA-II used for this work is its 
ability to model only a portion of the full cylindrical 
geometry for problems with cyclic symmetry. The grid used to 
model the slanted slot geometry comprised only a 30 degree 
sector of the full 360 degree geometry. This reduces the 
number of computational cells required in each case by a 
factor of 12 in comparison with a full 360 degree grid. 

The z and theta spacing of the k and j planes 
respectively are selected so that diagonal rows of vertices 
lie along the desired slot angle. The theta increment is 
defined as an angle and is referred to as thdiv below. The 
slot is divided up into 15 increments with respect to the j 
direction. 

L sina + W 
THDIV = R (n"y5 

m 

L - length of slot W g - slot width 

s 

a - slot angle R m - mixer radius 


With the theta increment selected the z increment is 
determined via trigonometry and is defined here as DZDIV. 


DZDI.V 


THDIV • R 
m 

tana 

s 
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With this type of mesh development the arc length of 
the modeled slot in the azimuthal direction corresponds 
exactly to the true 45 degree slot, but a slight discrepancy 
with the axial length of the actual and modeled slot exists. 
Vertex velocities are set to the jet velocity in groups of 
five in a k-plane and as the k-plane number is advanced 
( moving axially down the mixing zone ) the starting vertex 
of the inflow group is staggered with respect to the j 
direction as shown in figure (2-2) . The complete slanted 
slot encompasses 11 k-planes in the axial direction. Total 
axial length of the modelled slot is calculated below: 

Axial length = (11) (0.20259 cm) = 2.228 cm 
The true length of a 45 degree slot is: 

True length = (1.25 in) (2 . 54) (sin45) = 2.245 cm 
The percent error in our modeled axial length is: 

2.245 cm - 2.228 cm 

% difference (100) 

2.245 cm 

= 0.7570 % or less than 1 % 
Two types of velocity boundary conditions must be 
specified in KIVA-II. These are velocities specified at 
verticies and velocities specified at the 
cell-face centers. The latter are associated with the 
calculation of fluxing volumes used with the convective 
terms in the governing equations. In subroutine BCFC the dot 
product of left cell face area vector and face centered 
velocity is calculated (ual(i4)). Multiplying the ual(i4) 
quantity for cells lying along the slot edge by 2/3, since 
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these cells only have 3 of 4 vertices set with the inflow 
velocity, reduces the flux passing through the cell face to 
50% of that of cell faces lying fully within the slot. The 
following equations further describe the relationship 
described above: 

“f ■ ? ( “vl + S v2 + “v3 + 3 v4> 

At a cell face along the slot edge: 

S f = I (° + 3. + + Sj) - I u. 

For a cell face completely within the slot: 

5V = u^ • X^At 


For a cell face along the slot edge: 

6V = dr) (j- u . ) • it. At = \ u . • it. At 

i 4 J 1 2 3 1 

u f - face centered velocity. 

u v . - face vertex velocities (i=l,4). 

u. - jet velocity. 

it-^ - cell left face area vector. 

At - time step. 


6V - volume fluxed through the left 
cell face in a given cycle. 


The section of subroutine BCFC performing this function is 


listed below. 

v f 4 * 

a bs ( 

rt Inf l J f k ) ) 

v f 3 = 

a os ( 

rtlnf l J + l,k ) ) 

v f 7 = 

aos t 

rcl nf l J+l tk+1 ) ) 

yjf 3 = 

abs t 

rtlnf I J » k 1 j ) 


flxfac*vf3+vf4*vf7+vf3 

ualnxp*0«25*i (u( 1 3 ) *v f 3 I4)*vf4+u(l7)*vf7+u( 18). *vf8)*alx( U) 

♦ (v(l3)*vf3+v(i4)*vf4+v( I7)*vf7+v( l 8)*vfS)*alyl i4) 

f Ixf ac=cvigt t 2./3* » f lxf ac/4 . t flxfac.eq.3* ) 1 dIz<U)I 

uall Ul*cvijt(f lxfac*ualnxp, ual(U), i.aq.nxp .and. rtcalf.na.O 


QQM11220 
QUM112 JO 
QQMU240 
QQIU1250 
QQM11260 
QQHH270 
QQH11230 
UJM112V0 
QQHliiOO 
IQQH1131Q 
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The computing mesh for a 30 degree slot geometry is 
similar to that for the 45 degree slanted slot with a few 
exceptions. The method used to generate the slanted slot 
meshes yields optimum results with a slot at 45 degrees from 
the horizontal, which had an axial length error of less than 
1 percent. The number of divisions with respect to the 
azimuthal and axial direction determines the error between 
actual and modeled slot measurements. An error in either the 
axial or azimuthal direction of 8% was found to exist with 
most combinations. With this in mind the number of divisions 
within the 30 degree slot was picked to remain the same as 
the number used with the 45 degree slot. The 30 degree 
slot's THDIV and DZDIV components are calculated below. 


THDIV 


(1.25 in) (2 . 54) (sin30) + (0.3125 in) (2.54) 
(2.5 in) (2.54) (15 div) 


= 0.025 rad ( 180/3.14159 ) = 1.43239 degrees 


DZDIV 


(0.025rad) (2.5 in) (2.54) 
tan3 0 


= 0.27496 cm 

Modeled axial length = (11) (0 . 27496) = 3.0245 cm 
Actual Axial Length = (cos30) (1.25") (2.54) = 2.7496 cm 

3.0245cm — 2.7496cm 

% Difference (100) = 9.9% 

2 . 7496cm 


Another difference in grid geometry between the two 
cases investigated is that 5 additional cells in the 
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azimuthal direction were added to the 30° slot geometry. 

This was necessary since the slot itself encompassed only 
21.48' degrees and the sector mesh encompasses 30 degrees. 
This increase in cells yielded a mesh of 51,282 cells in 
comparison with the 39,627 cell, 45 degree baseline mesh. 

II. 4. Outlet Boundary Modifications 

The actual experimental rapid mixer test rig shown in 
figure (1-1) is equipped with a 6 ft lean zone. Since the 
flowfield becomes fully developed before the actual outflow 
port, a shorter computational lean zone length is used which 
also reduces execution run time. A 44.2 in lean zone was 
selected for the open end mesh. This long length ensured the 
bulk swirling profile dissipated prior to reaching the 
outflow boundary. 

Early in the transient of the nonreacting runs using 
the pressure outlet boundary conditions provided in the 
KIVA-II code flow reversal occurred on the outflow boundary. 
Once this flow reversal occurred, air at ambient 
temperature would enter the mesh. This phenomenon would 
greatly increase the required cpu time to reach steady 
state. Several options are available to prevent this flow 
reversal from occurring, two of which are presented with 
this work. By deactivating some of the exit plane cells, the 
outflow port is modified such that flow reversal is 
prevented. With a 18x16x86, 45 degree slanted slot mesh the 
first 6 cells in the i direction were left open and the 
remaining 12 were blocked out. Results from this approach 
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are presented to compare with a solution obtained by 
modifying the ouflow boundary conditions to prevent the flow 
reversal. This i-n effect leaves an open orifice in the 
center of a cylindrical blockage. This left a 2.857 cm open 
radius within the 8.89 cm radius lean zone. Figures (2-12) 
and (2-13) show the region of interest ( mixing section ) of 
the grid. This mesh used a shorter lean zone of only 
24.149 in than the open end mesh. History plots showing 
arrival of steady state with this partially blocked off mesh 
at J* = 20, are shown in figs. (2-3), (2-4), and (2-5). 
Figures (2-3) , (2-4) and (2-5) have cgs units. The partial 

blockage causes a 5% pressure increase throughout the mesh 
accompanied by a 5% decrease in the axial flow speed, when 
compared to results of a completely open end calculation. 
Results from the open end mesh are discussed in chapter 4 . 
Figures (2-6) , (2-7) , (2-8) , (2-9) , (2-10) , and (2-11) show the 
results of the partially blocked grid run. Qualitatively 
comparing the jet flow field characteristics with the open 
end mesh, the results compare well. A vortex occurring 
directly to the right, with respect to axial flow direction, 
(CCW) of the slanted slot can be seen in figure (2-10) . 
Secondly the low pressure region occurring directly 
downstream of the slot near the outer wall is also noted 
shown in figure (2-7) . A disadvantage of the partially 
blocked end grid is a greatly increased run time. The 
pressure increase as steady state is approached requires 
roughly 3 times more runtime than a similar run with an open 
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end grid. Average time steps from the J* = 20 blocked grid 
are compared with the J* = 25 open end grid below. 

(blocked) average dt = 0 . 1439s/9820cyc = 1.465e-05 s/cyc 
(open) average dt = 0 . 1067s/4693cyc = 2.273e-05 s/cyc 

The blocked mesh has a lower time step due to high shear 
stress at the partial blockage on the outflow boundary. 

To allow use of the open end grid while preventing any 
flow reversal along the outflow boundary required 
modifications to the outflow boundary conditions. Within 
subroutine BC a condition was set that would zero the 
velocity of any vertices with a negative value along the nzp 
plane, occurring only during the false transient. Also, in 
subroutine BCFC the face velocity would be set to zero along 
the nzp plane if the negative values occurred indicating 
flow reversal. Together these two modifications prevented 
flow reversal from occurring at any point during the 
transient. Steady state could be reached with 3000 to 4000 
time steps typically. 

II. 5. Turbulent Kinetic Energy Solution 

Another numerical anomaly experienced with KIVA-II 
during the nonreactive mixer runs was the formation of a 
localized tke high spot along the axis directly opposite the 
slanted slot. Knowing that with this geometry the turbulence 
should be at a maximum in the vicinity of the jet (slanted 
slot) [11] the localized high spot along the centerline was 
a numerical problem. The localized turbulent kinetic energy 
high spot is a result of insufficient accuracy with the 
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convergence criteria with the velocity and pressure 
solutions . 

The localized high spot would form almost immediately at 
the beginning of the transient and its magnitude was on the 
order of 100 times that of other surrounding cells. Once the 
localized high spot formed it would slowly migrate down the 
length of the mesh and not advect across the outflow 
boundary. The partially blocked end grid was able to pass 
the tke high spot out of the mesh. The quicker running 
simpler open end grid was unable to pass it through. The 
high spot trapped at the outflow would increase the 
amplitude of pressure oscillations and iterations for tke 
and turbulence dissipation. One run in which the localized 
high spot had formed at the outflow boundary and drove up 
the number of iterations for the turbulence solution was 
stopped at cycle 2749. Figures (2-14) , (2-15) , (2-16) , (2-17) , 
and (2-18) show plots from a run with high localized tke at 
the outlet centerline. Figure (2-14) of the velocity field 
looks well developed. The normal velocity plot in figure 
(2-15) also is well developed and can be compared with 
normal velocity plots from chapter 4 as well as the previous 
section of this chapter. Figures (2-17) and (2-18) reveal 
the localized tke high spot clearly. The tke high spot was 
always accompanied with a low pressure disturbance at the 
same location. Excerpts from file 12, and file 59 listing 
the iterations per cycle from a run with no tke anomaly 
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nearing steady state and the 2749 cycle with the turbulence 

problem are listed below. 

Localized TKE along the outlet- centerline. .. . 

ncyc= 2700 ca= 0.0 t= 1.3002E-01 dc= 4 . 3866E-05 gr=1.388 dcpuodt- 6.7358E+ 

0j c8h 1 8= 0 . 0000E+00 pm= 0.0000E+00 avp= 1.0136E+07 perr= 3.4771E-03 grow 

iters: v« 0 v= 8 t= 15 p= 10 k= 76 e= 73 big= 2 nfluxs= 7 pgs= 1.9442E 

+00 

urtin: -0.48100E+04 
raom-racio: 0.25009E+02 

ncyc = 2725 ca= 0.0 t= 1.3114E-01 dt= 4.4658E-05 gr=1.451 dcpuodt- 6.9149E+ 

0j c8hl8= 0 . OOOOE+OO pm- 0.0000E+O0 avp= 1.0130E+07 perr- 1 .4054E-03 rstr 

iters: y« 0 v- 13 t- 16 p- 7 k- 36 e- 50 big- 2 nfluxs- 7 pgs- 1.9453E 

400 

Normal rapid mixer run nearing steady state.... 

nc yc= 3299 ca- 0.0 t= 1.5447E-01 dt- 4.4655E-05 gr= 1.062 dcpuodt- 5.0615E+ 

0 c8hl8= 0 OOOOE+OO pm= 0. OOOOE+OO avp= 1.0155E+07 perr- 1.9643E-04 rstr 

iters: y- 0 v- 3 t* 7 p« 2 k= 6 e= 6 big= 2 nfluxs- 7 pgs- 1.9107E 

+00 

dtcon- 6 . 69053E-06 dcacc- 1.23886E-04 dtrsc- 4.46499E-05 dc- 4.46499E-05 
ncyc- 3300 ca= 0.0 t= 1.5451E-01 dt- 4.4650E-05 gr= 1.072 dcpuodt- 5.1069E+ 

0j c8hl8- 0. OOOOE+OO pm- 0. OOOOE+OO avp- 1.0155E+O7 perr- 1.3972E-04 rstr 

iters: y- 0 v- 3 t- 7 p- 2k- 6 e- 6 big- 2 nfluxs- 7 pgs- 1.9111E 

+00 

Initially, to eliminate the numerical anomaly a type of 
tke clipper was inserted into the code. It determined the 
average tke considering only the first 4 cells from the 
centerline along the entire length of the mesh. Then each 
cell within the 4 cell sweep range was compared with the 
average tke value and set to the average tke if exceeding 
it. Averaging was used because physically tke values near 
the centerline should be rather uniform. This technique 
prevented localized high spots from forming on the 
centerline with out altering the jet analysis within the 
area of interest ( k=l to k=60 ) . 
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Figure 2-la 
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BASELINE RAPID MIXER GEOMETRY 



Figure 2-lb 
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Figure 2-6 
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Figure 2-7 
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Figure 2-8 
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Figure 2-10 
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Figure 2-15 
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CHAPTER III 


PROCEDURE 

III.l. Initialization -of Run Conditions 

After comparing the characteristics and cpu 
requirements of the partially blocked mesh with the fully 
open mesh the latter was chosen for results generation. 
This mesh would take typically 3000 time steps from 
time = 0.0 sec to reach steady state. This corresponds to 
0.1-0.14 sec to reach the steady state. Table 3.1 lists 


several important parameters 
mixer runs. 

Table 


Variable 

Settina 

LWALL 

1 

JSECTR 

l 

IREZ 

0 

CYL 

1.0 

PGSSW 

1.0 

DTI 

1 . 04e-06 

THSECT 

30.0 

EPSY 

1.0e-03 

EPSV 

1.0e-03 

EPSP 

1.0e-04 

EPST 

1.0e-03 

EPSK 

1. 0e-03 

EPSE 

1.0e-03 

TCYLWL 

811.0 

PARDON 

0 

ANC4 

0.0 

ADIA 

0.0 

TKEI 

1600.0 

SGSL 

0.0 

Distamb 

0.0 


set m the input file for the 


Description 

selects law-of-the-wall 
boundary conditions 
sector mesh desired 
eulerian grid 
cylindrical mesh 
pressure gradient scaling 
method utilized 
initial time step (s) 
sector mesh of 30 degrees 
convergence mass diffusion 
convergence momentum diff. 
convergence of pressure 
convergence of heat diff. 
convergence for tke diff. 
convergence of epsi diff. 
wall temp = 811° kelvin 
quasi second order upwind 
differencing 
no node coupling 
isothermal walls 
2.5% of inlet velocity 
use k- e turb. model 
specified outlet pressure 
directly at nzp plane 
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III. 2. Inf low /Outflow Boundaries 

The bottom boundary is entirely open for inflow, while 
the right boundary allows inflow only along the slot and is 
a solid curved surface elsewhere. Air at 1000° F enters the 
computing mesh across the bottom boundary, while air at 
77° F enters through the slanted slot. Pressure and 
temperature of the entering air are defined by specifying 
their stagnation state and using the following insentropic 
relations : 


SPDIN(M) = SPDINO (M) ( — 


SPDINO(M) - the stagnation density for species 
m set in the input file. 

P - pressure of the respective cell into 

which the species is being fluxed. 

gamma - ratio of the specific heats determined 

at the respective inflow boundary conditions. 


PO 


- stagnation pressure at the respective 
inflow boundary set in the input file. 


Along the bottom boundary turbulent length scale is set 
to 0.065, and the tke itself is based upon w' (fluctuating 
velocity) specified as 2.4% of the square of the mean flow 
velocity. Stagnation pressure and temperature are set to 10 
atm and 1000° F, respectively. 

Along the right inflow boundary the stagnation pressure 
and temperature are set to 10 atm and 77° F, respectively. 
Three percent of the perpendicular distance across the slot 
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opening (0.2213 in) is used to calculate both tke and length 
scale entering with the jet. 

A static pressure is applied along the outlet boundary 
by setting DISTAMB to 0.0 in the input file. DISTAMB is a 
quantity that allows the applied pressure boundary to be 
moved away from the physical outflow boundary, which in some 
cases reduces time to reach steady state. The pressure 
applied is 10 atm. Acoustic wave reflection was found not to 
be a problem. It was seen that giving DISTAMB a positive 
value increased the probability of flow reversal at the 
outlet boundary. 

III. 3. Initial Problem Conditions 

Species densities of the individual components of air 
at 10 atm and 1000° F were calculated and set in the input 
file to initialize the respective densities through out the 
mesh at time = 0.0 sec. The initial velocity profile was set 
based upon a relation with the square of the radius at a 
given axial position in the geometry. This approach reduced 
required time to reach steady state. The reference velocity 
of Win and radius of the rich zone are used in the following 
relation: 

W(i4) = Win*(- fffi ) 2 

W(i4) - value of axial velocity at time=0.0 sec. 

set at every vertex during initial setup. 

Win - the axial velocity applied to the bottom 
inflow boundary at k=l. 

R(k) - radius at a given k-plane. 
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This type of initial velocity distribution provided good 
results in comparison with runs using an initial velocity of 
0.0 through out the mesh. Another alternative is ramping the 
axial inflow velocity ( win ) linearly, but this required a 
longer transient than runs using the described initial 
distribution . 

III. 4. Jet Development 

The jet was ramped from an initial velocity of 0.0 to 
the desired value of the variable URTIN in the input file. A 
major advantage of setting the jet velocity in the input 
file is that subsequent higher values could be set and the 
run resumed from a previous J* steady state binary file via 
KIVA-II's restart capability. This decreased transient time 
for the new higher J* steady state by roughly a third from 
the transient time required for a run starting at time = 0.0 
seconds. 

The jet is ramped very slowly early in the transient to 
minimize pressure waves between the right boundary and the 
centerline. Between cycle 1 and 120 the jet is increased by 
0.01 cm/s per le-05 sec. From cycle 120 to 800 the ramp is 
increased to 1.0 cm/s per le-05 sec. Ramping the jet based 
upon time and not cycle early in the transient is 
advantageous because the time step (dt) fluctuates and drops 
into the low (1.0e-06 sec. . . 1 . 2e-06 sec) range at some 
points. Beyond cycle 800 the jet ramp is raised to 5 cm/s 
per cycle and up to 10 cm/s per cycle after cycle 1000. 
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III. 5. Steady State Determination 

Steady state was determined from visual inspection of 
history plots produced after a restart. This file contains 9 
different histories (3 pressure, 3 temperature, 3 velocity) 
from three different locations in the mesh. Appendix A lists 
the specific history locations for all calculations 
including the partially blocked mesh described in chapter 2. 
It was observed that steady state would first be established 
first at the k = 1 plane and last at the outlet. Since the 
important region of the grid is in the vicinity of the jet, 
histories indicating steady state at the end of the region 
of interest (k=60) were used to determine acceptable 
development of the jet and termination of the run. When one 
of the history locations was changed to the mesh outlet 
(k=nz) for the 30 degree slanted slot runs, it was clearly 
evident that temperature was not steady at cycle 3000 and 
required several thousand more cycles to approach steady 
state. Allowing the entire mesh to reach steady state was 
not realistic with every run due to excessive cpu time 
requirements, thus only the region of interest was 
considered for steady state. Appendix B contains history 
plots at the end of the region of interest used to verify 
the flow transient had reached steady state for the result 
cases presented in chapters 4 and 5. 
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CHAPTER IV 


45° SLANTED SLOT RAPID MIXER RESULTS 
JV.l. Initial Setup 

The baseline mixer utilizes a 45 deg slanted slot 
orifice to inject 77° F secondary air into the mixing 
section of the rapid mixer. The slot has a 4 to 1 aspect 
ratio. The results presented in this chapter were produced 
with a 20x16x110 open ended mesh utilizing the tke averaging 
technique discussed in chapter 2. Modifications to the 
outflow boundary conditions of the KIVA-II code set any 
components of negative face or vertex velocities to zero. 
This is only necessary during the transient and does not 
effect the steady state results, since there is no tendency 
for the flow field to reverse after the jet has been ramped 
to is full inflow speed. The mesh is shown in full length in 
figure (4-1) . Figure (4-2) shows the azimuthal cell spacing 
in the j direction. Figures (4-3) and (4-4) are enlargements 
of the region of interest (k=l to k=60) . This mesh models 
44.2 in of the lean zone which has an actual length 
of 6 ft. 

IV. 2. Bulk Swirl Profile 

An important characteristic of the flowfield, which 
enhances thermal mixing, is the swirl profile originating at 
the 45 degree slanted slot which persists well into the lean 
zone. Figures (4-5) to (4-11) are k-slice velocity field 
plots with J*=5, that show the swirling velocity profile 
along the length of the mixer geometry. Figures 
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(4-12) to (4-18) are the temperature contour plots at the 
same axial locations. 

Air injected from the slanted slot remains along the 
outer perimeter surrounding a high temperature inner core. 

At k=4 2 the low temperature contour line is distorted from 
effects of jet penetration radially inwards into the mixer. 
This is a full 5.27 cm from the downstream edge of the 
slanted slot. At k=55, which is 17.57 cm from the downstream 
edge of the slanted slot, the low temperature contour is 
smoothed from turbulent mixing between the high and low 
temperature air masses. At k=65 and beyond, the temperature 
contours assume a circular profile. At higher J* values the 
high temperature contour line will intersect the centerline. 

The jet interacts with the axial mainstream inducing a 
CCW spin in the region near the wall. An opposite CW spin 
develops between this outer CCW swirling profile and the 
centerline. Figure (4-6) in the 21° expansion region of the 
grid reveals that as the axial velocity decreases in the 
expansion the magnitude of swirl velocity increases. A well 
developed swirl profile is seen at k=55, but at k=65 the 
swirling profile is beginning to vanish. At k=75,85 and 95 
the azimuthal component of velocity has decreased greatly 
with the entire profile assuming a CW swirl of a small 
magnitude. Figures (4-19) and (4-42) of normal velocity 
contours reveal that the defined swirl profile extends only 
slightly beyond the area of interest (k=60) . 
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IV. 3. Higher Momentum Flux Ratios 

As jet velocity increases, the jet penetrates further 
into the mixing region. 'The vortex directly to the right of 
the jet shown in figures (4-21) through (4-23) elongates 
considerably. Temperature and turbulent kinetic energy 
contours shown in figures (4-24) through (4-29) reveal 
increased radial penetration of the jet. 

At low momentum flux ratios a stagnant low pressure 
region exists just downstream of the jet near the wall seen 
in figure (4-33) . As J* increases the low pressure point 
moves away from the wall following the penetration of the 
jet (figure (4-34)). In Fig. (4-35) the low pressure region 
has shifted to the center of the mixer geometry, which 
indicates the jet is in effect pinching the axial flow. 

The corresponding temperature contours are shown in figures 
(4-36) through (4-38) . The high temperature contour line 
moves almost directly radially inward intersecting the 
centerline in Fig. (4-38) . Figures (4-39) through (4-41) 
reveal high turbulence near the wall at the base of the jet 
as expected in this region of high shear. 

Increased momentum flux ratio increases the extent of 
the swirling flow as shown in figures (4-45) through (4-50) , 
which correspond to two axial locations. 
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IV. 4. Numerical Accuracy 

The mesh used for the results presented in the previous 
section consists of 39,627 cells. This allows for adequate 
spatial resolution to resolve the high shear region in the 
vicinity of the jet. A run with J =25 was conducted with a 
coarser mesh to verify grid independence of the generated 
results. The coarser cell mesh used 16,660 cells (13x16x69) 
and is shown in figure (4-63) . Figures (4-51) to (4-56) 
generated with the coarse mesh compare well with figures 
(4-57) to (4-62) . Temperature contour values generated with 
the coarser mesh agree within 4%. Comparing the velocity 
plots a 4% difference in u and v velocity components, and an 
11% difference in w velocity components is noted. The 
pressure contour plots show agreement within 8%. All 
quantities compare well with the exception of the turbulent 
kinetic energy, which is 40-60% lower in the coarser mesh. 
This is most likely a result of the tke averaging employed 
with both runs. The tke is averaged from i=l to i=4 along 
the centerline in both runs. The coarse grid having 7 fewer 
cells in the radial region has a greater portion of the same 
geometry averaged. At J =120 the temperature contour in 
Fig. (4-26) does not indicate any coalescence of the jets. 
This is a characteristic of quasi-second-order-upwind 
differencing. When using partial and full donor cell 
differencing coalescence is observed at much lower momentum 
flux ratios. 
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JV.5. Lean Zone Length Effects 

A J*=60 run was conducted with a mesh having 20 less 
k-planes and an axial length decrease- of 27.94 cm from the 
20x16x110 grid. Results from the run are presented with the 
equivalent plots from the 44.2 in lean zone grid in figures 
(4-64) to (4-71) . Visual qualitative comparison of plots 
from the region of interest of both grids reveals very 
little deviation. Variation of flowfield parameters shown in 
these figures is tabulated in table 4.1. 

TABLE 4.1 


Plot 


Velocity 


Variation 


UM 

0.53% 

L-Contour 


VM WM 

3.59% 1.21% 

H-Contour 


Normal Velocity 2.84% 3.59% 

Pressure 0.87% 0.078% 

Temperature 0.05% 0.27% 

TV. 6. Bulk Temperature Profile 

Bulk temperature is one characteristic that is of 

interest when considering nonreactive flow mixing 

performance of the rapid mixer geometry. Bulk temperature is 

defined by the following equation: 

/pc uTdA 

Bulk Temperature = — 

/pc udA 
v 

p - fluid density u - fluid velocity 

T - temperature A - area 

c v - specific heat at constant volume 
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Bulk temperature is the temperature to which a non- 
isothermal fluid mixture would reach if allowed to come to 
an equilibrium state. Bulk temperature profiles for the 
20x16x110 grid at J=5,25,40,60 and 120 are presented in 
figure (4-72) . The bulk temperature is plotted verses mixer 
length within the area of interest. The runs were 
terminated prior to the entire lean zone reaching steady 
state due to high cpu time requirements. 

Bulk temperature was calculated on each k-plane of the 
sector mesh once the area of interest reached steady state. 
The bulk temperature plummets at the point of jet injection. 
The minimum bulk temperature is lower with increasing 
momentum flux ratio. 

TV. 7. Swirl Number Profile 

Beer and Chigier [12] proposed the non-dimensional 
parameter Swirl Number. This parameter characterizes the 
amount of rotation in a given axial flow. Swirl number is 
defined by the following equations: 


G = /2nr$rpl5dr 
m 

G t = /2nrp3^dr + /2nrpdr 
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tangetial velocity 
radius 

Swirl number = G /G. 

m t 

This quantity like bulk temperature was calculated on each 
k-plane once the area of interest had reached steady state. 

Lefebvre [1] states most fixed vane swirlers of 
practical interest operate in the range of (SN>0.6). The 
fixed vane swirler physically surrounds the fuel injector 
and its purpose is to induce rotation of the secondary air 
which creates a recirculation region within the primary 
zone . 

The swirl number profiles for all five momentum flux 
ratios are plotted verses axial location in figure 
(4-73) . Except for a spike at the point of jet injection the 
swirl number remains between 3.0e-04 and -5.0e-04, which is 
essentially zero compared to 0.6 in conventional gas 
turbines utilizing fixed vane swirlers. The swirl profile 
indicates the conservation of angular momentum within the 
flowfield. Angular momentum of the outer CCW rotating field 
is virtually offset by the inner CW rotating core. 


U - axial velocity W - 

p - static pressure r - 

p - fluid density 
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Figure 4-34 
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Figure 4-42 
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CHAPTER V 


30° SLANTED SLOT RESULTS 
V.l. 30° Slanted Slot Geometry 

As discussed in chapter 2 the 30° slanted slot was 
modeled in a similar fashion to the 45° slot. This approach 
required an additional 5 j -planes, which increased the mesh 
size to 51,282 cells. Figures (5—1) and (5—2) show the mesh 
within the region of interest ( k=l to k=60 ) , and figure 
(5-3) shows azimuthal spacing of the j— planes. Effect of 
momentum flux ratio variation was conducted using 
J =5,25,40,60 and 120 at a fixed density ratio of 2.76. 

These runs also utilized tke averaging along the mesh 
centerline from i=l to 4 to prevent the development of 
localized turbulent kinetic energy high spots. 

V.2. Jet Description 

The slanted slot jet interaction with the axial 
mainstream is similar to that for the 45° slanted slot. 
Holdeman et al. [10] theorized 45° slots produce less 
penetration and mixing performance than equivalent area 
parallel slots. The 30° slanted slot having less slant 
produces greater radial penetration into the mixing region 
of the mesh. Jet penetration can be visually characterized 
in figures (5—4) to (5—9) . Comparing these velocity and 
temperature profiles across the slanted slot with 
corresponding figures in chapter 4 the increased penetration 
is readily evident. The 30° slanted slot geometry having 
increased jet penetration develops a low pressure region 
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just downstream of the jet along the centerline and at lower 
momentum flux ratios than the 45° slot geometry. The low 
pressure region is developed at J*=40 while it was not 
evident until J*=60 with the 45° slot. Regardless of the 
difference in jet penetration, the bulk temperature 
profiles shown in figure (5-23) drop to values similar to 
the profiles generated with the 45° slot. 

V.3. Swirl Profile 

The 30° slots having a smaller inclination with the 
centerline, deflects less of the mainstream into the outer 
CCW swirl region. Accordingly the inner CW rotating bulk 
swirl profile encompasses a greater portion of the cross 
sectional area. 

The j -plane normal velocity plots shown in figures 
(5-19) to (5-21) clearly indicate a reduced magnitude of the 
azimuthal velocity in the outer CCW profile in comparison 
with the 45° slanted slot. An increase in magnitude of 
azimuthal velocity in the inner CW profile is most evident 
at higher momentum flux ratios. 

As with the 45° slanted slot, the swirl number 
variation with axial distance in fig. (5-22) , reveals a net 
cancelation of angular momentum. These values in comparison 
with those referenced by Lefebvre [1] are effectively zero. 
V.4. tke Averaging with the 30 ° Slot 

The averaging out of high tke spots near the centerline 
is one solution to remedy this problem. This method works 
well if the jet does not penetrate into the region 
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undergoing tke averaging ( i=l to i=4 ) . When the jet 
physically enters this region its interaction with the 
mainstream results in high turbulence that tke averaging 
alters. This is a problem with the 30° slanted slot at 
J*=120. A simpler method preventing the numerical high tke 
high spots is decreasing velocity convergence criteria to 
1.0e-04 and pressure convergence criteria to l.0e-05. This 
alternative eliminates the need to "smooth” the centerline 
tke high spots and is recommended in general. 
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Figure 5-2 
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Figure 5-3 
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Figure 5-7 
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Figure 5-11 
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Figure 5-14 
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Figure 5-15 
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Figure 5-19 
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Figure 5-20 
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Figure 5-21 
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CHAPTER VI 


CONCLUSIONS & RECOMMENDATIONS 
VI. 1. Conclusions and Flowfield Evaluation 

Results presented in chapter 4 and 5 were produced 
utilizing a tke averaging technique from i=i-4 along the 
centerline. This method yields the best results with low to 
mid momentum flux ratios, because the tke values are nearly 
uniform around the centerline. At high momentum flux ratios 
(J =120) when the jet penetrates near the centerline, this 
assumption is no longer valid. By superimposing values of 
turbulence lower than physically realistic, turbulent 
diffusion is decreased and an unrealistic central channel 
with excessively high magnitudes of velocity develops along 
the centerline. This region can be seen on the velocity 
profiles for the J* = 120 runs shown in chapter 4 and 5. it 
is also manifested in the temperature and tke contours. 
Overall development and trends within the flowfield are 
characterized below: 


1 . 


Jet radial penetration 
momentum flux ratios. 


increases with higher 


2. The bulk swirling profile extends deeper 
into the lean zone with higher J* values. 

3 . Jet radial penetration increases with 
decreasing slot angle. 

4. Net angular momentum from the 
inner CW and outer CCW profiles is 
virtually zero as shown by the 
dimensionless parameter swirl number. 

This is consistent with the conservation 
of angular momentum. 
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VI .2 . 


Future Work 


Small ripples in the bulk temperature profiles from 
chapter 4 -and 5 are related to tke averaging. Bulk 
temperature profiles generated with runs not using tke 
averaging along the centerline were found to be much 
smoother. The disadvantages of tke averaging are more 
pronounced when investigating temperature difference between 
high and low contour values as a function of axial location 
down stream of the point of jet injection. Significant jumps 
in temperature difference were noted from one axial location 
to the next. It is clear the tke averaging along the 
centerline is unacceptable if accurate evaluation of thermal 
mixing performance is to be conducted. 

To eliminate the need to average tke, convergence 
criteria must be decreased to l.0e-05 for the pressure 
solution and to 1.0e-04 for the velocity solution. Runs 
conducted at lower J* values show a time step 2 to 3 times 
higher than the corresponding run utilizing tke averaging. 
Figures (6-1) to (6-5) present flowfield results from a 
baseline steady state run with J* = 120 and the reduced 
convergence criteria. These figures are much improved from 
the corresponding plots shown in chapters 4 and 5. Two 
important questions are yet to be answered from this work: 

1 * =^narai sianted siot promote better mixing than 
a parallel slot. 
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2. Can the experimental results presented 

in a URTC report, which state that optimum 
mixing with^the baseline geometry occurs at 
mid-range J values, be validated with the 
computational model. 

Holdeman [10] has concluded that slots with a 45° slant 
in a rectangular duct produce reduced penetration and poorer 
thermal mixing than parallel slots. Fox and Gex [13] 
conducted an experimental investigation of the effect of the 
angle of injection of a single jet into a large reservoir of 
fluid, while attempting to minimize induced swirl and 
vorticity. Oldshue [14] states swirl is detrimental to jet 
mixing. It is unclear at this point if swirl within the 
circular axisymmetric baseline geometry is detrimental to 
mixing performance. If the jet and hot axial mainstream have 
mixed perfectly at a given z coordinate downstream of the 
point of jet injection there will be no variation of 
temperature along the cross section. The profile will be 
flat and have the numerical value equivalent to bulk 
temperature at the given z coordinate. Future work will 
include calculations for various momentum flux values with 
45° slot and parallel slot geometries, and a comparison of 
temperature departure from the bulk mean temperature as a 
function of axial position. These calculations will take 
advantage of improved convergence criteria for velocity and 
pressure solutions to eliminate the need for tke averaging 
along the centerline. 
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APPENDIX A 


CASE GEOMETRY REVIEW 


Definition of Geometric Variables : 

ZM - Z coordinate of the point of transition 

between the 30° contraction and the mixing 
region. 

ZBOC - Z coordinate of the upstream edge of the 
slanted slot. 

ZEOC - Z coordinate of the downstream edge of the 
slanted slot in. 

ZML - Z. coordinate of the transition between the 
mixing region and the 21° expansion. 

ZLN - Z coordinate of the transition between the 
21° expansion and the lean zone. 


45 Degree Slanted Slot 39,627 Cell Mesh: 


ZM = 
ZBOC = 
Z(k=25) = 
ZEOC = 
Z (k=40) = 

ZML = 
ZLN = 
Z (k=55) = 
Z (k=60) = 
Z (k=65) = 
Z (k=75) = 
Z (k=85) = 
Z (k=95) = 


7.28 cm 
13.15 cm 
14.03 cm 
15.38 cm 

19.44 cm 
21.25 cm 
27.87 cm 
32.95 cm 
39.30 cm 
45.65 cm 
58.34 cm 
74.05 cm 

99.45 cm 


j = 9 Corresponds to 14.624° 


Lean Zone Length = 44.18" included with this mesh. 
Corresponding Figures: 


2-1,2-2,2-14,2-15,2-16,2-17,2-18, 
4-1 to 4-50,4-57 to 4-62,4-65,4-67, 
4-69,4-71,4-72,4-73 
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30° Slanted Slot 51,282 Cell Mesh: 

ZM = 7.179 cm 
ZBOC = 12.88 cm 
Z(k=25) = 13.99 cm 
ZEOC = 15.64 cm 
Z(k=40) = 20.62 cm 
ZML = 21.25 cm 
ZLN = 27.87 cm 
Z(k=55) = 35.49 cm 
Z(k=60) = 41.84 cm 

j = 9 Corresponds to 11.459° 

Lean Zone Length = 45.183" included with this mesh. 
Corresponding Figures: 

5-1 to 5-23 


45° Slanted Slot Partially Blocked Off 28,101 Cell Mesh: 

ZM = 7.28 cm 
ZBOC = 13.15 cm 
Z (k=25) = 14.37 cm 
ZEOC = 15.38 cm 
ZML = 21.25 cm 
ZLN = 27.87 cm 
Z(k=60) = 39.04 cm 

j = 9 Corresponds to 14.624° 

Lean Zone Length = 24.15" included with this mesh. 
Corresponding Figures: 

2-3 to 2-13 
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45 0 Slanted Slot 16,660 Cell Mesh: 


ZM 

ZBOC 

Z(k=28) 

ZEOC 

ZML 

ZLN 

Z(k=60) 


7.28 cm 
13.15 cm 

14.87 cm 
15.38 cm 
21.25 cm 

27.87 cm 
37.02 cm 


j = 9 Corresponds to 14.624° 


Lean Zone Length = 29.99" included with this mesh. 
Corresponding Figures: 


4-51 to 4-56,4-63 


45° Slanted Slot Reduced Length 32,487 Cell Mesh: 
Lean Zone Length = 33.20" included with this mesh. 


a. Results generated with tke averaging along the 
centerline. 

Corresponding Figures: 

4-64,4-66,4-68,4-70 

b. Results generated without tke averaging utilizing 
increased convergence criteria for the velocity 
and pressure solution. 

Corresponding Figures: 

6-1 to 6-5 
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APPENDIX B 


History Plots 


Steady state determination was facilitated with three 
sets of history plots for each geometry. Each set included 
pressure, temperature and velocity. The most upstream 
location is found to reach steady state first and the 
furthest downstream location last. If pressure, temperature 
and velocity are at steady state along the nz plane ( the 
mixer outlet ) it can be assumed the flowfield has reached 
steady state at all other regions of the mesh. To conserve 
cpu required to reach steady state the entire lean zone was 
not allowed to reach steady state prior to terminating the 
run. Instead, history plots from a point slightly downstream 
the transition point between the 21° expansion and the lean 
zone were used to determine if the area of interest ( the 
mixing region ) had reached steady state. The following are 
i,j and k coordinates where the history plots were taken: 

45° Slanted Slot Mesh: i,j,k = 8,8,60 
45° Slanted Slot Coarse Mesh: i,j,k = 7,8,60 
30° Slanted Slot Mesh: i,j,k = 11,8,60 

These locations correspond to positions at the end of the 
area of interest and roughly centered with respect to the 
radial and azimuthal directions within the 30 degree sector 
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mesh. Table B.l categorizes history plots included with this 
appendix by mesh type. 


TABLE B.l 

45° Slanted Slot Mesh 

B-l , B-2 , B-3 
B-4 , B-5 , B-6 
B-7 , B-8 , B-9 
B-10 , B-ll , B-l 2 
B-13 , B-14 , B-15 


30° Slanted Slot Mesh 

J* = 5,25,60,120 B-16 , B-17 , B-18 

J = 40 B-l 9 , B-20 , B-21 


45° Slanted Slot Coarse Mesh 
J * = 25 B-2 2 , B-2 3 , B-2 4 



J* = 60 
J = 120 


167 



168 







Figure B-2 
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Figure B-3 
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Figure B-4 
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Figure B-7 
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Figure B-8 
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Figure B-9 
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Figure B-10 
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Figure B-14 
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Figure B-17 
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Figure B-19 
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2, Computational M ethod 

A gas turbine combustor is a complex 
system in which a number of physical 
processes occur simultaneously. These 
processes include subsonic fluid flow, 
turbulent transport of mass, momentum, 
and energy, wall heat transfer, turbulent 
flame propagation, formation and quench- 
ing, dilution jet mixing, pollutant for- 
mation, and perhaps fuel spray dynamics. 
The overall combustor performance is the 
net result of a complex interplay among 
these various effects. Extensive appli- 
cation of the ICEd-ALE finite volume 
method, primarily to internal combustion 
engines, have demonstrated it to be a 
reliable and robust technique for the 
modeling of chemically reactive fluid 
dynamics with fuel sprays. ^”7 

The ICEd-ALE technique is a hybrid 
of the Implicit Cont inuous- f lu id Eule- 
rian (ICE) 8 and Arbitrary Lagrangian- 
Eulerian ( ALE methods . The former pro- 
vides flow-speed versatility through an 
implicit treatment of the pressure cal- 
culation while the latter allows for 
arbitrarily shaped hexahedronal shaped 
cells which may be used to model curved 
or moving boundaries. The overall cal- 
culation is a result of contributions 
from separate Eulerian and Lagrangian 
phases. This tends to result in reduced 
numerical diffusion in contrast to purely 
Eulerian difference methods. 

An overview of the governing equa- 
tions solved by the ICEd-ALE finite 
volume scheme follows. 

Conservation Equations 

Mass : 

^ + ?.(p3) = + p l 

o t in mm 

for the mth species, where is the mass 
density of^ species m, u is the fluid 
velocity, J m is the ^i f f us ive mass flux 
of species m, and p m is the rate of 
change of P m due to chemical reaction. 
The diffusive mass flux is given by 
Fick’s Law: 


3 = pD^Y 

m m 


where P is the total mass density and 
is the species mass fraction. The mass 
diffusivity D, which is assumed to be the 
same for all species, is determined using 
the Schmidt number Sc, (assumed 
cons tan t ) : 


D 


P 

pSc 


where p is viscosity. 


Momentum : 

- + 3.(pSS) ■ -^p - ^pk 4* ^ • q 

where p is the fluid pressure, and k is 
the specific turbulent kinetic energy. 
The viscous stress tensor is Newtonian 
and is expressed as: 

<2 = y[?u + ($S) T ] + X$-uU 


where P and ^ are the shear and bulk 
viscosity coefficients, U is the unit 
dyadic tensor, and superscript T denotes 
the transpose. The viscosity coeffi- 
cients are given by: 


P “ 


air 


4 - c p 
P 



e 


and 



where e is the specific turbulent kinetic 
energy dissipation rate and c is an 
empirical turbulence constant equalling 
0.09. A Sutherland formula is used for 
the molecular viscosity of air with A^ 
and A^ constant: 

. 3/2 


air 


V 

T + 


*2 


Energy : 

d 4 * (pul) ■= - $-q 4 -pe 4 -Q c 


where I is the specific internal energy 
and Q is a source term due to chemical 
heat release. The total heat flux q, is 
gained by summing the effects of heat 
conduction and enthalpy transport by mass 
diffusion : 

q “ q + q. = -k$t - zh 3 

^ d m m m 

where K is the thermal conductivity, T is 
the fluid temperature, and is the 
species specific enthalpy. The thermal 
conductivity is determined using the 
Prandtl number Pr (assumed constant): 



where Cp is the specific heat at constant 
pressure . 

Turbulence Equations 

A two-equation k-e turbulence model 
is used. It solves the following pair of 
equations for k and z : 

pk^*u = 


^ + *-< P Sk) + 1 


+ $*(p 5 ) = o 


o:^u 4- Vk] - pe 


2 



TF + + f| c + c ) 3 = 


■c pe. 

L 2 


Ipr ^ C 1 + f l c r o : $u-< 
c K E ! ' 

where c C j, c e2 , c c 3 , |‘r k , and !t £ are 

vaTue r s^lO C ° nStan, - s the following 


The global reaction rate expression is: 
• P b ’ o a 

“ ' • ■ = k fr m (or) - k t _ n( ») 


c = 1 . < 1-1 
1 


= 1 . 9 2 


C r = - 1.0 

♦ 3 


K, 

br 


Pr k = 1.0 Pr £ = 1.3 

Jenith^caYe 1 °" “ """ E 1 the turbulent 
length scale L is computed directly from: 

k 3 ' 2 


L = c 


u e e 


where 


1/2 


p c= r 

FTFT 7 1 • 

S tate and Chemica l EouaHnn. 

. ., ThC s t flte relations are assumed to 
be, those of an ideal gas mixture, there- 

P " V fiCp./w.). 

I(T) = m ^ P m^ p ) (T) 

C P (T) = £ (P m /p ) C pm< T >. 

h m (T > = V T > + R 0 T/W o • 


thermo? iS , the univer sal gas constant- w 
the molecular weight of species me ^ 

s^creTmf 10 at press^e^ 

«- r - 

is used to model th- u C s mechanism 
fuel-air mixture: * COmbustlon of the 

2 C 3 H 8 + 70 2 + 8 H 2 ° + 6C0 
2C0 + 0 2 •<-+ 2C0 

0 + N 2 ++ N + NO 
0 + NO +-*■ o + n 
N + OH H + NO 


global "frm V W 7 “ *^ br “ 1 {J- 

whe re a- «„ d b* are no” necessarilv 1 1, 

sto.chiometric coefficients b». 

r b ;rr 

are determined from the Arrhen^s: ‘ 

k = ATCe xp(fPf) 

The temperature exponent £, the active 
lion energy E n anH tv, active- 

term A i ’ ? d the Pre-exponen t i al 

tion with 1S ° determin ed from correla- 
with experimental results. 

Six equilibrium reactions are considered: 

H 2 -<-* 2H 

0 2 20 

N 2 -*->• 2N 

0 2 + H 2 -<-*• 20H 

°2 + 2 H 2 ° A0H 

0 2 + 2C0 2C0 2 . 


The rate of each equilibri 
-rmined implicitly from: 


dete 


um reaction is 


S <r> 


K c (T) 


coef f icienta^of V spe a ^ e ^ stoichi °metric 
right side, resp^ctivei; i 0 , L n ?% le ( f T t I " d 
the known CQuilibrinm c is 

- * ^ uor „ v b t r .“" »*>->■ 

— Computational Model 


Computing 

equal [y r spaced abound' 1 th^rap^n 

peSfr :-- 

direction. Hence, a 30 deg 6 
dimensional sector ».ck g ‘ three- 

320 cells is used fo^ . < ?° mposed °f 32, 

mesh consists of 20 cells m V Tl ' ^ 

vertices) in , . 7 1 planes of 

cells (17 i — , i radla l (r) direction, 16 
us U7 j -planes of vertices) 
azimuthal f 0) din^^t- r cices; i n the 
(102 ni erection, and 101 cells 

llu * k-planes of vertices) in + u ls 

Is) direction IFig 2 , t,. , ' “ xi * 1 

iurrzx. r sh "Mr 
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Boundary Conditions 

The mainstream inlet {3) conditions 
and jet inlet (J) conditions are both 
assumed to be uniform. Since the j = 1 
and j = 17 azimuthal vertex planes are 
periodic the radial and azimuthal velo- 
city components must be equal at corres- 
ponding vertices on these planes. Simi- 
larly, convective fluxing is adjusted to 
ensure equality between efflux and influx 
of any scalar variable on corresponding 
cell surfaces. Specified temperatures 
are set for and T for mainstream and 
jet inlets. Turbulent law-of- the-wall 
formulations are used to represent the 
velocity and thermal boundary layers at 
solid walls. The inflow turbulent kine- 
tic energy is set as 10% of the mean 
inflow kinetic energy: 

k 3 = .l(V*/2) kj = . l(Vj/2) 

where is the mainstream inlet axial 
velocity and Vj is the jet radial inlet 
velocity. The inflow turbulent length 
scales are : 


.1 R 


,1 W 


3 "r “J s 

where R r is the rich zone radius and W 
is the slot width. s 

The slot inflow was specified on 
square cell surfaces on the i = 21 plane 
in the slot region ( k = 35 - 46). The 
slot slant was represented by specifying 
V j along diagonal rows of vertices (Fig. 
3). For cells bisected by the slot edge 
the mass flow rate was reduced by one 
half of that for cells lying completely 
within the slot. 

4 ♦ Results 
Nonreactive Flow 

A comparison is made between slots 
parallel to the combustor centerline with 
slots slanted at 45 deg. The flow con- 
ditions are : 

P 3 ” 10 atm = 1.013 x lO^dynes/cm^ f 
T 3 = 811 k, V 3 - 1600 cm/s, 

Tj = 298 k 

The jet inlet velocity Vj was determined 
for a given j e t- t o -ma i ns t ream momentum 
flux ratio : 

p j v j 

j* = J J 

D V 2 
M 3 3 

Results were determined for values of J* 

= 4 0 and 6 0 . 

For the parallel slots, the effect 
of an axial crossflow on a purely radial 
injection jet produces a pair of counter- 


rotating vortices with mirror symmetry 
between the slots when viewed on a given 
z-plane (Fig. 4), For slanted slots, 
however only a single vortex develops 
between slots (Fig. 5). For slots 

slanted in a clockwise sense, ail the 
vortices will have a clockwise rotation. 
These twelve vortices near the wall of 
the quick-quench mixer act in effect as 
swirlers spaced around the circumference 
of the mixer. The cumulative effect is 
to produce a net clockwise swirl velocity 
incipient at the tailing slot edge and 
persisting downstream into the lean zone 
region. The outer clockwise swirl velo- 
city induces an inner counterclockwise 
swirl velocity region between the center- 
line axis and the outer swirl region 
(Fig. 6) The downstream swirl develop- 
ment is shown in Fig. 7. Note that the 
outer swirl region is clearly distin- 
guished by positive contour lines which 
terminate at the combustor wall, while 
the inner swirl region has negative con- 
tour lines which tend to form a closed 
region. This type of swirl velocity 
pattern has been observed experimentally 
for slanted slot jet injections into a 
cylindrical crossflow. 

The relative mixing effectiveness 
between the parallel slot and slanted 
slot cases was evaluated in terms of 
deviation from the bulk internal energy 
at a given axial (z) location. The bulk 
internal energy passing through a given 
cross section is: 


where m is the total mass flow rate 
through the cross section. The mixing 
effectiveness (n) is then defined as: 


/ .1.1-1, l<ta 


Z b“ 

The variation of n with z is shown in 
Fig. 8 for both the slanted slots (indi- 
cated by (45)) and parallel slot injec- 
tion . 

Reactive Flow 

One reactive flow results case is 
presented which corresponds to equiva- 
lence ratios of 1.6 and .5 in the rich 
and lean zones respectively. The opera- 
ting conditions are: 

P 3 = 6.89 x 10^ dynes/ cm 2 (100 psia) , 

m f = 63 • 9 g/s. m ar = 625.9 g/s, 

T, = 811 k 
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where m f and m ar are the fuel and rich 
zone air mass flow rates respectively. 

The slanted slots once again produce 
a swirl velocity development downstream 
of the slot region (Figs. 9-10), The jet 
in effect acts as a flameholder. The 
consumption of propane and fuel mass 
fraction are shown in Figs. 11-12. Tem- 
perature contours are shown in Figs. 13- 
1*1. The momentum flux ratio based on the 
rich zone momentum flux is J* r 28. The 
average cross sectional temperature along 
the combustor is shown in Fig. 15, The 
pattern factor along the combustor is 
shown in Fig. 16. The pattern factor is 
de f i ned as : 

T - f 
PF = max 4 

wh t T < : T3 

where T 3 is the average inlet tempera- 
ture, is the average temperature of a 

cross section, and l m ax is the maximum 
temperature of a cross section. Emission 
index variation for CO and NO is shown in 
Figs. 17-18. 


5 . Cone 1 us 1 ons 

The IbEd-ALE method is a versatile 
and accurate computational method for the 
modeling of the reactive flowfields 
inside gas turbine combustors. Radial 
dilution air jet injection through 
slanted slots into a mainstream flow in a 
duct of circular cross sectional area 
results in a bulk swirl velocity 
downstream of the slot region. The swirl 
velocity considerably enhances thermal 
mixing in comparison to parallel slot 
injection. A completely preraixed fuel/air 
supply was used to minimize thermal 
N0 X production in the rich zone and 
isolate the N0 X contribution from the 
quick- quench mixer and lean zone regions. 
Excessive NO x generation < exit emission 
index 1.2 ) was not observed for the case 

presented. Further analysis, along with 
experimental validation, is necessary for 
a more thorough unders tand ing of the RQL 
combustion processes. 
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Figure 2. Computational mesh. 
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Figure 7. Swirl velocity contours, H=5116 cm/s, L=-2517 cm/s (j=9,J*=60). 



Figure 9. Swirl velocity contours, H=12277 cm/s, L=-3963 cm/s (j=9,J*=28). 
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Figure 10. Swirl velocity (k=60 , J*=28) . 


Figure 11. Propane depletion. 



Figure 12. 


Propane mass fraction, H=.0286, 


L-. 00318, (j=9, J*=28) . 
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Figure 17. NO emission index 



emission index 


Figure 16. Pattern factor. 


Figure 18. CO 


